rm(list=ls())

options(stringsAsFactors = FALSE)

library(Hmisc)
library(ggplot2)
library(corrplot)
library(car)
library(stargazer)
library(reshape2)

# > sessionInfo()
# R version 4.2.0 (2022-04-22)
# Platform: aarch64-apple-darwin20 (64-bit)
# Running under: macOS Monterey 12.4
# 
# Matrix products: default
# LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
# 
# locale:
#   [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] reshape2_1.4.4  stargazer_5.2.3 car_3.0-13      carData_3.0-5   corrplot_0.92   Hmisc_4.7-0    
# [7] ggplot2_3.3.6   Formula_1.2-4   survival_3.3-1  lattice_0.20-45

load("data_psys.Rda")
load("data_psys_splits.Rda")

# Table 2. Descriptive statistics ----

data_for_desc <- data_psys[, c("closure_mean", "alt_mean", "form_mean", "acc_mean", "v2x_libdem_mean", "v2x_polyarchy_mean", "v2x_liberal_mean", "v2xcl_rol_mean", "v2x_jucon_mean", "v2xlg_legcon_mean", "enpp_mean", "pol_mean", "tev_mean", "tor_mean", "log10_gdp_mean", "age")]

stargazer(data_for_desc,
          type = "text",
          digits = 2)

# Figure 1. Correlations of the indices and sub-indices of closure and liberal democracy in Europe. ---- 

corrplot_data <- data_psys[, c("closure_mean", "alt_mean", "form_mean", "acc_mean", "v2x_libdem_mean", "v2x_polyarchy_mean", "v2x_liberal_mean", "v2xcl_rol_mean", "v2x_jucon_mean", "v2xlg_legcon_mean")]

names(corrplot_data) <- c("closure", "alternation", "formula", "access", "liberal dem.", "electoral dem.", "liberal comp.", "rule of law", "judicial", "legislative")

corr_out <- rcorr(as.matrix(corrplot_data),type="pearson")

corrplot(corr_out$r, method="color", col="white",  
         type="upper",
         addCoef.col = "black", # Add coefficient of correlation
         tl.col="black", #Text label color and rotation
         # Combine with significance
         p.mat = corr_out$P, sig.level = 0.05, 
         # hide correlation coefficient on the principal diagonal
         diag=FALSE,
         cl.pos = FALSE,
         insig = "blank",
         addgrid.col = "black"
)

# Table 3. Relationship between closure and measures of democracy. ----

m00 <- lm(v2x_libdem_mean ~ closure_mean +  enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m0 <- lm(v2x_polyarchy_mean ~ closure_mean +  enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m1 <- lm(v2x_liberal_mean ~ closure_mean +  enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m2 <- lm(v2xcl_rol_mean ~ closure_mean +  enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m3 <- lm(v2x_jucon_mean ~ closure_mean +  enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m4 <- lm(v2xlg_legcon_mean ~ closure_mean +  enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)

m00_std_data <- data.frame(scale(m00$model))
m0_std_data <- data.frame(scale(m0$model))
m1_std_data <- data.frame(scale(m1$model))
m2_std_data <- data.frame(scale(m2$model))
m3_std_data <- data.frame(scale(m3$model))
m4_std_data <- data.frame(scale(m4$model))

m00 <- lm(v2x_libdem_mean ~ closure_mean +  enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m00_std_data)
m0 <- lm(v2x_polyarchy_mean ~ closure_mean +  enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m0_std_data)
m1 <- lm(v2x_liberal_mean ~ closure_mean +  enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m1_std_data)
m2 <- lm(v2xcl_rol_mean ~ closure_mean +  enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m2_std_data)
m3 <- lm(v2x_jucon_mean ~ closure_mean +  enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m3_std_data)
m4 <- lm(v2xlg_legcon_mean ~ closure_mean +  enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m4_std_data)

models <- list(m00,m0,m1,m2,m3,m4)

stargazer(models,
          type="text",
          #report=('vc*p'),
          omit.stat=c("f", "ser"),
          star.cutoffs = c(0.05, 0.01, 0.001)
)

# Table 4. Relationship between the components of closure and measures of democracy. ----

m00 <- lm(v2x_libdem_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m0 <- lm(v2x_polyarchy_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m1 <- lm(v2x_liberal_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m2 <- lm(v2xcl_rol_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m3 <- lm(v2x_jucon_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m4 <- lm(v2xlg_legcon_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)

m00_std_data <- data.frame(scale(m00$model))
m0_std_data <- data.frame(scale(m0$model))
m1_std_data <- data.frame(scale(m1$model))
m2_std_data <- data.frame(scale(m2$model))
m3_std_data <- data.frame(scale(m3$model))
m4_std_data <- data.frame(scale(m4$model))

m00 <- lm(v2x_libdem_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m00_std_data)
m0 <- lm(v2x_polyarchy_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m0_std_data)
m1 <- lm(v2x_liberal_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m1_std_data)
m2 <- lm(v2xcl_rol_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m2_std_data)
m3 <- lm(v2x_jucon_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m3_std_data)
m4 <- lm(v2xlg_legcon_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m4_std_data)

models <- list(m00,m0,m1,m2,m3,m4)

stargazer(models,
          type="text",
          omit.stat=c("f", "ser"),
          star.cutoffs = c(0.05, 0.01, 0.001)
)

# Variance inflation factors

lapply(models, car::vif)

# Outliers and influential cases

influencePlot(m00)
influencePlot(m0)
influencePlot(m1)
influencePlot(m2)
influencePlot(m3)
influencePlot(m4)

# Appendix. Table 1. Relationship between the components of closure and measures of democracy, interactions with region. ----

m00 <- lm(v2x_libdem_mean ~ region_ew + alt_mean*region_ew + form_mean*region_ew + acc_mean*region_ew + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m0 <- lm(v2x_polyarchy_mean ~ region_ew + alt_mean*region_ew + form_mean*region_ew + acc_mean*region_ew + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m1 <- lm(v2x_liberal_mean ~ region_ew + alt_mean*region_ew + form_mean*region_ew + acc_mean*region_ew + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m2 <- lm(v2xcl_rol_mean ~ region_ew + alt_mean*region_ew + form_mean*region_ew + acc_mean*region_ew + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m3 <- lm(v2x_jucon_mean ~ region_ew + alt_mean*region_ew + form_mean*region_ew + acc_mean*region_ew + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m4 <- lm(v2xlg_legcon_mean ~ region_ew + alt_mean*region_ew + form_mean*region_ew + acc_mean*region_ew + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)

m00_std_data <- data.frame(scale(m00$model))
m0_std_data <- data.frame(scale(m0$model))
m1_std_data <- data.frame(scale(m1$model))
m2_std_data <- data.frame(scale(m2$model))
m3_std_data <- data.frame(scale(m3$model))
m4_std_data <- data.frame(scale(m4$model))

m00 <- lm(v2x_libdem_mean ~ region_ew + alt_mean*region_ew + form_mean*region_ew + acc_mean*region_ew + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m00_std_data)
m0 <- lm(v2x_polyarchy_mean ~ region_ew + alt_mean*region_ew + form_mean*region_ew + acc_mean*region_ew + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m0_std_data)
m1 <- lm(v2x_liberal_mean ~ region_ew + alt_mean*region_ew + form_mean*region_ew + acc_mean*region_ew + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m1_std_data)
m2 <- lm(v2xcl_rol_mean ~ region_ew + alt_mean*region_ew + form_mean*region_ew + acc_mean*region_ew + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m2_std_data)
m3 <- lm(v2x_jucon_mean ~ region_ew + alt_mean*region_ew + form_mean*region_ew + acc_mean*region_ew + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m3_std_data)
m4 <- lm(v2xlg_legcon_mean ~ region_ew + alt_mean*region_ew + form_mean*region_ew + acc_mean*region_ew + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m4_std_data)

models <- list(m00,m0,m1,m2,m3,m4)

stargazer(models,
          type="text",
          omit.stat=c("f", "ser"),
          star.cutoffs = c(0.05, 0.01, 0.001)
)

# Appendix. Table 2. Relationship between the components of closure and measures of democracy, with split party systems (n=67). ----

m00 <- lm(v2x_libdem_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys_splits)
m0 <- lm(v2x_polyarchy_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys_splits)

m1 <- lm(v2x_liberal_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys_splits)
m2 <- lm(v2xcl_rol_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys_splits)
m3 <- lm(v2x_jucon_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys_splits)
m4 <- lm(v2xlg_legcon_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys_splits)

m00_std_data <- data.frame(scale(m00$model))
m0_std_data <- data.frame(scale(m0$model))
m1_std_data <- data.frame(scale(m1$model))
m2_std_data <- data.frame(scale(m2$model))
m3_std_data <- data.frame(scale(m3$model))
m4_std_data <- data.frame(scale(m4$model))

m00 <- lm(v2x_libdem_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m00_std_data)
m0 <- lm(v2x_polyarchy_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m0_std_data)
m1 <- lm(v2x_liberal_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m1_std_data)
m2 <- lm(v2xcl_rol_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m2_std_data)
m3 <- lm(v2x_jucon_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m3_std_data)
m4 <- lm(v2xlg_legcon_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m4_std_data)

models <- list(m00,m0,m1,m2,m3,m4)

stargazer(models,
          type="text",
          omit.stat=c("f", "ser"),
          star.cutoffs = c(0.05, 0.01, 0.001)
)

# Appendix. Table 3. Relationship between the components of closure and measures of democracy, using ministerial volatility instead of alternation.  ----

m00 <- lm(v2x_libdem_mean ~ mv_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m0 <- lm(v2x_polyarchy_mean ~ mv_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m1 <- lm(v2x_liberal_mean ~ mv_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m2 <- lm(v2xcl_rol_mean ~ mv_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m3 <- lm(v2x_jucon_mean ~ mv_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m4 <- lm(v2xlg_legcon_mean ~ mv_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)

m00_std_data <- data.frame(scale(m00$model))
m0_std_data <- data.frame(scale(m0$model))
m1_std_data <- data.frame(scale(m1$model))
m2_std_data <- data.frame(scale(m2$model))
m3_std_data <- data.frame(scale(m3$model))
m4_std_data <- data.frame(scale(m4$model))

m00 <- lm(v2x_libdem_mean ~ mv_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m00_std_data)
m0 <- lm(v2x_polyarchy_mean ~ mv_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m0_std_data)
m1 <- lm(v2x_liberal_mean ~ mv_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m1_std_data)
m2 <- lm(v2xcl_rol_mean ~ mv_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m2_std_data)
m3 <- lm(v2x_jucon_mean ~ mv_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m3_std_data)
m4 <- lm(v2xlg_legcon_mean ~ mv_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m4_std_data)

models <- list(m00,m0,m1,m2,m3,m4)

stargazer(models,
          type="text",
          omit.stat=c("f", "ser"),
          star.cutoffs = c(0.05, 0.01, 0.001)
)

# Appendix. Table 4. Relationship between the components of closure and measures of democracy, with proportions of no and wholesale alternation. 

m00 <- lm(v2x_libdem_mean ~ none + wholesale + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m0 <- lm(v2x_polyarchy_mean ~ none + wholesale + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m1 <- lm(v2x_liberal_mean ~ none + wholesale + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m2 <- lm(v2xcl_rol_mean ~ none + wholesale + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m3 <- lm(v2x_jucon_mean ~ none + wholesale + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)
m4 <- lm(v2xlg_legcon_mean ~ none + wholesale + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys)

m00_std_data <- data.frame(scale(m00$model))
m0_std_data <- data.frame(scale(m0$model))
m1_std_data <- data.frame(scale(m1$model))
m2_std_data <- data.frame(scale(m2$model))
m3_std_data <- data.frame(scale(m3$model))
m4_std_data <- data.frame(scale(m4$model))

m00 <- lm(v2x_libdem_mean ~ none + wholesale + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m00_std_data)
m0 <- lm(v2x_polyarchy_mean ~ none + wholesale + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m0_std_data)
m1 <- lm(v2x_liberal_mean ~ none + wholesale + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m1_std_data)
m2 <- lm(v2xcl_rol_mean ~ none + wholesale + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m2_std_data)
m3 <- lm(v2x_jucon_mean ~ none + wholesale + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m3_std_data)
m4 <- lm(v2xlg_legcon_mean ~ none + wholesale + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m4_std_data)

models <- list(m00, m0, m1,m2,m3,m4)

stargazer(models,
          type="text",
          omit.stat=c("f", "ser"),
          star.cutoffs = c(0.05, 0.01, 0.001)
)

# Appendix. Table 5. Relationship between the components of closure and measures of democracy, influential cases excluded (n=49).

toinclude <- c("Albania","Austria I","Austria II","Belgium","Bulgaria", "Croatia","Cyprus", "Czechia",  "Denmark",  "Estonia I","Estonia II","Finland I","Finland II", "France II","France III","France IV","Georgia","Germany I", "Germany II", "Greece II","Greece IV","Hungary",  "Iceland","Ireland","Italy", "Latvia I", "Latvia II","Lithuania","Luxembourg", "Moldova", "Netherlands","North Macedonia", "Norway", "Poland I", "Poland II","Portugal I","Portugal II","Romania", "Serbia", "Slovakia", "Slovenia", "Spain I","Spain II", "Spain III","Sweden", "Switzerland","Turkey I", "Turkey II", "United Kingdom" )

data_psys_sub <- data_psys[data_psys$country %in% toinclude, ]

m00 <- lm(v2x_libdem_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys_sub)
m0 <- lm(v2x_polyarchy_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys_sub)
m1 <- lm(v2x_liberal_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys_sub)
m2 <- lm(v2xcl_rol_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys_sub)
m3 <- lm(v2x_jucon_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys_sub)
m4 <- lm(v2xlg_legcon_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = data_psys_sub)

m00_std_data <- data.frame(scale(m00$model))
m0_std_data <- data.frame(scale(m0$model))
m1_std_data <- data.frame(scale(m1$model))
m2_std_data <- data.frame(scale(m2$model))
m3_std_data <- data.frame(scale(m3$model))
m4_std_data <- data.frame(scale(m4$model))

m00 <- lm(v2x_libdem_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m00_std_data)
m0 <- lm(v2x_polyarchy_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m0_std_data)
m1 <- lm(v2x_liberal_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m1_std_data)
m2 <- lm(v2xcl_rol_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m2_std_data)
m3 <- lm(v2x_jucon_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m3_std_data)
m4 <- lm(v2xlg_legcon_mean ~ alt_mean + form_mean + acc_mean + enpp_mean + pol_mean + tev_mean + tor_mean + log10_gdp_mean + age, data = m4_std_data)

models <- list(m00,m0,m1,m2,m3,m4)

stargazer(models,
          type="text",
          omit.stat=c("f", "ser"),
          star.cutoffs = c(0.05, 0.01, 0.001)
)

# Appendix. Figure 1. Sensitivity to the exclusion of single cases. ----

countries <- unique(data_psys$country)

DVs <- c("v2x_libdem_mean", "v2x_polyarchy_mean", "v2x_liberal_mean", "v2xcl_rol_mean", "v2x_jucon_mean", "v2xlg_legcon_mean")

jack_models <- list()

for (i in 1:length(countries)) {
  country <- countries[i]
  data_jack <- data_psys[data_psys$country != country,]
  models_out <- list()
  for (j in 1:length(DVs)) {
    formul <- formula(paste(DVs[j], "~", paste(c("alt_mean", "form_mean", "acc_mean", "enpp_mean", "pol_mean", "tev_mean", "tor_mean", "log10_gdp_mean", "age"), collapse = "+"), sep = ""))
    out <- lm(formul, data = data_jack)
    std_data <- data.frame(scale(out$model))
    out <- lm(formul, data = std_data)
    models_out[[j]] <- out
  }
  jack_models[[i]] <- models_out
}

names(jack_models) <- as.character(countries)

test <- lapply(jack_models, function(x) lapply(x, function(y) summary(y)))
test_pval <- lapply(test, function(x) lapply(x, function(y) coef(y)[,3]))

test_pval <- lapply(test_pval, function(x) do.call(rbind, x))

test_pval <- as.data.frame(do.call(rbind, test_pval))

test_pval$Country <- rep(countries, each = length(DVs))

test_pval$DV <- rep(DVs, length(countries))

test_plot <- melt(test_pval, id = c("Country", "DV"))

test_plot <- test_plot[test_plot$variable != '(Intercept)', ]

test_plot$DV <- factor(test_plot$DV, levels = DVs, labels = c('DV: liberal democracy', 'DV: polyarchy', 'DV: liberal component', 'DV: rule of law', 'DV: judicial constraints', 'DV: legislative constraints'))

test_plot$variable <- car::recode(test_plot$variable, "'alt_mean' = 'alternation';'form_mean' = 'formula';'acc_mean' = 'access';'enpp_mean' = 'fragmentation';'pol_mean' = 'polarization';'tev_mean' = 'volatility';'tor_mean' = 'type of regime';'log10_gdp_mean' = 'log10 GDP';'age' = 'age'")

test_plot$variable <- factor(test_plot$variable, levels = c('alternation', 'formula', 'access', 'fragmentation','polarization','volatility','type of regime','log10 GDP','age'))

ggplot(test_plot, aes(x = variable, y = value)) +
  #geom_point() +
  geom_text(aes(label = Country), alpha=0.5, cex = 1.5) +
  geom_hline(aes(yintercept = 2), color = "red", lty = "dashed") +
  geom_hline(aes(yintercept = -2), color = "red", lty = "dashed") +
  geom_hline(aes(yintercept = 1.67), color = "blue", lty = "dotted") +
  geom_hline(aes(yintercept = -1.67), color = "blue", lty = "dotted") +
  facet_wrap(~DV) +
  labs(title = NULL,
       subtitle  = 'The figure shows the t-values of model coefficients when a particular country is excluded from the analysis. Dashed lines\nindicate the thresholds for statistical significance at the 0.05 level and dotted lines significance at the 0.1 level.') +
  ylab("t-value") +
  xlab(NULL) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Appendix. Figure 2. Sensitivity to the exclusion of single cases, influential cases excluded (n=49). ----

countries <- unique(data_psys_sub$country)

DVs <- c("v2x_libdem_mean", "v2x_polyarchy_mean", "v2x_liberal_mean", "v2xcl_rol_mean", "v2x_jucon_mean", "v2xlg_legcon_mean")

jack_models <- list()

for (i in 1:length(countries)) {
  country <- countries[i]
  data_jack <- data_psys_sub[data_psys_sub$country != country,]
  models_out <- list()
  for (j in 1:length(DVs)) {
    formul <- formula(paste(DVs[j], "~", paste(c("alt_mean", "form_mean", "acc_mean", "enpp_mean", "pol_mean", "tev_mean", "tor_mean", "log10_gdp_mean", "age"), collapse = "+"), sep = ""))
    out <- lm(formul, data = data_jack)
    std_data <- data.frame(scale(out$model))
    out <- lm(formul, data = std_data)
    models_out[[j]] <- out
  }
  jack_models[[i]] <- models_out
}


names(jack_models) <- as.character(countries)

test <- lapply(jack_models, function(x) lapply(x, function(y) summary(y)))
test_pval <- lapply(test, function(x) lapply(x, function(y) coef(y)[,3]))

test_pval <- lapply(test_pval, function(x) do.call(rbind, x))

test_pval <- as.data.frame(do.call(rbind, test_pval))

test_pval$Country <- rep(countries, each = length(DVs))

test_pval$DV <- rep(DVs, length(countries))

test_plot <- melt(test_pval, id = c("Country", "DV"))

test_plot <- test_plot[test_plot$variable != '(Intercept)', ]

test_plot$DV <- factor(test_plot$DV, levels = DVs, labels = c('DV: liberal democracy', 'DV: polyarchy', 'DV: liberal component', 'DV: rule of law', 'DV: judicial constraints', 'DV: legislative constraints'))

test_plot$variable <- car::recode(test_plot$variable, "'alt_mean' = 'alternation';'form_mean' = 'formula';'acc_mean' = 'access';'enpp_mean' = 'fragmentation';'pol_mean' = 'polarization';'tev_mean' = 'volatility';'tor_mean' = 'type of regime';'log10_gdp_mean' = 'log10 GDP';'age' = 'age'")

test_plot$variable <- factor(test_plot$variable, levels = c('alternation', 'formula', 'access', 'fragmentation','polarization','volatility','type of regime','log10 GDP','age'))

ggplot(test_plot, aes(x = variable, y = value)) +
  #geom_point() +
  geom_text(aes(label = Country), alpha=0.5, cex = 1.5) +
  geom_hline(aes(yintercept = 2), color = "red", lty = "dashed") +
  geom_hline(aes(yintercept = -2), color = "red", lty = "dashed") +
  geom_hline(aes(yintercept = 1.67), color = "blue", lty = "dotted") +
  geom_hline(aes(yintercept = -1.67), color = "blue", lty = "dotted") +
  facet_wrap(~DV) +
  labs(title = NULL,
       subtitle  = 'The figure shows the t-values of model coefficients when a particular country is excluded from the analysis. Dashed lines\nindicate the thresholds for statistical significance at the 0.05 level and dotted lines significance at the 0.1 level.') +
  ylab("t-value") +
  xlab(NULL) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
